home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / vol_100 / 167_01 / spline.c < prev    next >
Text File  |  1985-11-13  |  18KB  |  735 lines
  1. /* 
  2. HEADER:     CUG
  3. TITLE:        SPLINE.C - Interpolate Smooth Curve;
  4. VERSION:    3.00;
  5. DATE:        09/26/86;
  6. DESCRIPTION:    "SPLINE takes pairs of numbers from the standard
  7.         input as abscissae and ordinates of a function.
  8.         It produces a similar set, which is approximately
  9.         equally spaced and includes the input set, on the
  10.         standard output. The cubic spline output has two
  11.         continuous derivatives and sufficiently many
  12.         points to look smooth when plotted.";
  13. KEYWORDS:    spline, graphics, plot, filter, UNIX;
  14. SYSTEM:        ANY;
  15. FILENAME:    SPLINE.DOC;
  16. WARNINGS:    NONE;
  17. CRC:        xxxx
  18. SEE-ALSO:    SPLINE.C;
  19. AUTHORS:    Ian Ashdown - byHeart Software;
  20. COMPILERS:    ANY;
  21. REFERENCES:    AUTHORS: Bell Laboratories;
  22.         TITLE:     UNIX Programmer's Manual, Volume 1
  23.              p. 145, Holt, Rinehart and Winston;
  24.         AUTHORS: R.W. Hamming;
  25.         TITLE:     Numerical Methods for Scientists and
  26.              Engineers, 2nd Edition
  27.              pp. 349 - 355, McGraw-Hill (1973);
  28.         AUTHORS: Forsythe, Malcolm & Moler;
  29.         TITLE:     Computer Methods for Mathematical
  30.              Computations
  31.              pp. 70 - 83, Prentice-Hall;
  32. ENDREF
  33. */
  34.  
  35. /*-------------------------------------------------------------*/
  36.  
  37. /* SPLINE.C - Interpolate Smooth Curve
  38.  *
  39.  * Version 3.00        September 26th, 1986
  40.  *
  41.  * Modifications:
  42.  *
  43.  *   V1.00 (85/11/01)    - beta test release
  44.  *   V2.00 (85/12/25)    - general revision
  45.  *   V2.01 (86/09/13)    - corrected command-line option bug for
  46.  *              Microsoft C compiler
  47.  *   V3.00 (86/09/26)    - added non-uniform abscissa spacing
  48.  *              capability 
  49.  *
  50.  * Copyright 1985,1986:    Ian Ashdown
  51.  *            byHeart Software
  52.  *            620 Ballantree Road
  53.  *            West Vancouver, B.C.
  54.  *            Canada V7S 1W3
  55.  *
  56.  * This program may be copied for personal, non-commercial use
  57.  * only, provided that the above copyright notice is included in
  58.  * all copies of the source code. Copying for any other use
  59.  * without previously obtaining the written permission of the
  60.  * author is prohibited.
  61.  *
  62.  * Synopsis:    SPLINE [option] ...
  63.  *
  64.  * Description:    SPLINE takes pairs of numbers from the standard
  65.  *        input as abscissae and ordinates of a function.
  66.  *        (A minimum of four pairs is required.) It
  67.  *        produces a similar set, which is approximately
  68.  *        equally spaced and includes the input set, on the
  69.  *        standard output. The cubic spline output (R.W.
  70.  *        Hamming, "Numerical Methods for Scientists and
  71.  *        Engineers", 2nd ed. 349ff) has two continuous
  72.  *        derivatives and sufficiently many points to look
  73.  *        smooth when plotted.
  74.  *
  75.  *        The following options are recognized, each as a
  76.  *        separate argument:
  77.  *
  78.  *        -a  Supply abscissae automatically (they are
  79.  *            missing from the input); spacing is given by
  80.  *            the next argument or is assumed to be 1 if
  81.  *            next argument is not a number.
  82.  *
  83.  *        -k  The constant "k" is used in the boundary
  84.  *            value computation
  85.  *
  86.  *            y " = ky " , y " = ky "
  87.  *             0    1     n         n-1
  88.  *
  89.  *            is set by the next argument. By default,
  90.  *            k = 0. A value of k = 0.5 often results in a
  91.  *            smoother curve at the endpoints than the
  92.  *            default value. Negative values for k are not
  93.  *            allowed. Cannot be used with -p option.
  94.  *
  95.  *        -n  Next argument (which must be an integer)
  96.  *            specifies the number of intervals that are to
  97.  *            occur between the lower and upper "x" limits.
  98.  *            If -n option is not given, default spacing is
  99.  *            100 intervals.
  100.  *
  101.  *        -p  Make output periodic, i.e. match derivatives
  102.  *            at ends. First and last input values must
  103.  *            agree. Cannot be used with -k option.
  104.  *
  105.  *        -x  Next 1 (or 2) arguments are lower (and upper)
  106.  *            "x" limits. Normally these limits are
  107.  *            calculated from the data. Automatic abscissae
  108.  *            start at lower limit (default 0). If either
  109.  *            argument is outside of the range of
  110.  *            abscissae, it is ignored.
  111.  *
  112.  * Diagnostics:    When data is not strictly monotone in "x", SPLINE
  113.  *        reproduces the input without interpolating extra
  114.  *        points.
  115.  *
  116.  * Bugs:    A limit of 1000 input points is silently
  117.  *        enforced.
  118.  *
  119.  *        The -n option has not been implemented in
  120.  *        accordance with the "UNIX Programmer's Manual"
  121.  *        specification. This was done to avoid ambiguities
  122.  *        when the -n option follows the -x option with one
  123.  *        argument.
  124.  *
  125.  *        At certain negative values for the -k option (for
  126.  *        example, k equals -4.0), the curve becomes
  127.  *        discontinuous. The -k option value has thus been
  128.  *        arbitrarily constrained to be greater than or
  129.  *        equal to zero.
  130.  *
  131.  * Credits:    The above description is a reworded and expanded
  132.  *        version of that appearing in the "UNIX Programmer's
  133.  *        Manual", copyright 1979, 1983 Bell Laboratories.
  134.  */
  135.  
  136. /*** Definitions ***/
  137.  
  138. #define FALSE         0
  139. #define TRUE         1
  140. #define MAX_SIZE  1000    /* Input point array limit */
  141.  
  142. #define ILL_ARG         0    /* Error codes */
  143. #define ILL_CMB      1
  144. #define ILL_KVL         2
  145. #define ILL_NVL         3
  146. #define ILL_OPT      4
  147. #define ILL_XVL         5
  148. #define INS_INP      6
  149. #define MIS_KVL      7
  150. #define MIS_NVL      8
  151. #define MIS_XVL      9
  152. #define MIS_YVL        10
  153. #define NMT_ORD        11
  154.  
  155. #define SQUARE(a) a*a
  156. #define CUBE(a) a*a*a
  157.  
  158. /*** Typedefs ***/
  159.  
  160. typedef int BOOL;    /* Boolean flag */
  161.  
  162. /*** Include Files ***/
  163.  
  164. #include <stdio.h>
  165. #include <ctype.h>
  166. #include <math.h>
  167.  
  168. /*** Main Body of Program ***/
  169.  
  170. int main(argc,argv)
  171. int argc;
  172. char **argv;
  173. {
  174.   int n = 0,
  175.       i,
  176.       j,
  177.       n_val = 0,
  178.       atoi();
  179.   float    x[MAX_SIZE],
  180.     y[MAX_SIZE],
  181.     h[MAX_SIZE-1];
  182.   double a_val = 1.0,
  183.      k_val = 0.0,
  184.      x1_val = 0.0,
  185.      x2_val = 0.0,
  186.      x_intvl,
  187.      ix,
  188.      iy,
  189.      d2y[MAX_SIZE],
  190.      atof(),
  191.      spl_int();
  192.   char buffer[257],
  193.        *temp,
  194.        *gets();
  195.   BOOL aflag = FALSE,    /* Command-line option flags */
  196.        kflag = FALSE,
  197.        pflag = FALSE,
  198.        x1flag = FALSE,
  199.        x2flag = FALSE,
  200.        is_float();
  201.   void spl_coeff(),
  202.        pspl_coeff(),
  203.        error();
  204.  
  205.   /* Parse the command line for user-selected options */
  206.  
  207.   while(--argc)
  208.   {
  209.     temp = *++argv;
  210.     if(*temp != '-')    /* Check for legal option flag */
  211.       error(ILL_OPT,*argv);
  212.     else
  213.     {
  214.       temp++;
  215.       switch(toupper(*temp))
  216.       {
  217.     case 'A':    /* "-a" option */
  218.       aflag = TRUE;
  219.       if(argc > 1 && is_float(*(argv+1)))
  220.       {
  221.         argc--;
  222.         argv++;
  223.          if((a_val = atof(*argv)) <= 0.00)
  224.           error(ILL_ARG,*argv);
  225.       }
  226.       break;
  227.     case 'K':    /* "-k" option */
  228.       if(pflag == TRUE)
  229.         error(ILL_CMB,NULL);
  230.       kflag = TRUE;
  231.       if(argc > 1 && is_float(*(argv+1)))
  232.       {
  233.         argc--;
  234.         argv++;
  235.         k_val = atof(*argv);
  236.         if(k_val < 0.00)
  237.           error(ILL_KVL,*argv);
  238.         break;
  239.       }
  240.       else
  241.         error(MIS_KVL,NULL);
  242.     case 'N':    /* "-n" option */
  243.       if(argc > 1)
  244.       {
  245.         argc--;
  246.         argv++;
  247.         if((n_val = atoi(*argv)) < 1)
  248.           error(ILL_NVL,*argv);
  249.         else
  250.           break;
  251.       }
  252.       else
  253.         error(MIS_NVL,NULL);
  254.     case 'P':    /* "-p" option */
  255.       if(kflag == TRUE)
  256.         error(ILL_CMB,NULL);
  257.       pflag = TRUE;
  258.       break;
  259.     case 'X':    /* "-x" option */
  260.       x1flag = TRUE;
  261.       if(argc > 1 && is_float(*(argv+1)))
  262.       {
  263.         argc--;
  264.         argv++;
  265.         x1_val = atof(*argv);
  266.       }
  267.       else
  268.         error(MIS_XVL,NULL);
  269.       if(argc > 1 && is_float(*(argv+1)))
  270.       {
  271.         x2flag = TRUE;
  272.         argc--;
  273.         argv++;
  274.         x2_val = atof(*argv);
  275.         if(x2_val <= x1_val)
  276.           error(ILL_XVL,x2_val);
  277.       }
  278.       break;
  279.     default:    /* "-n" option */
  280.       error(ILL_OPT,*argv);
  281.       }
  282.     }
  283.   }
  284.   if(n_val == 0)    /* Set "n_val" if not given */
  285.     n_val = 100;
  286.  
  287.   /* Get the input data */
  288.  
  289.   while(1)        /* ... while there is more input data */
  290.   {
  291.     if(aflag == TRUE)    /* Automatic abscissae were called for */
  292.     {
  293.       if(n == 0)
  294.     x[0] = x1_val;
  295.       else
  296.     x[n] = x[n-1] + a_val;
  297.     }
  298.     else        /* Abscissae supplied with input data */
  299.     {
  300.       if(gets(buffer))
  301.     x[n] = atof(buffer);
  302.       else
  303.     break;
  304.     }
  305.     if(gets(buffer))    /* Read in the corresponding ordinate */
  306.       y[n] = atof(buffer);
  307.     else
  308.     {
  309.       if(aflag == TRUE)
  310.     break;
  311.       else
  312.     error(MIS_YVL,NULL);
  313.     }
  314.     if(++n == MAX_SIZE)    /* Max